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SUMMARY 


Two new version of the k — uj two-equation turbulence model will be presented. The new Baseline 
(BSL) model is designed to give results similar to those of the original k-u> model of Wilcox, but without 
its strong dependency on arbitrary freestream values. The BSL model is identical to the Wilcox model in 
the inner 50% of the boundary-layer but changes gradually to the high Reynolds number Jones-Launder 
k-e model (in a k-u formulation) towards the boundary-layer edge. The new model is also virtually 
identical to the Jones-Launder model for free shear layers. The second version of the model is called Shear- 
Stress Transport (SST) model. It is based on the BSL model, but has the additional ability to account 
for the transport of the principal shear stress in adverse pressure gradient boundary-layers. The model is 
based on Bradshaw’s assumption that the principal shear-stress is proportional the turbulent kinetic energy, 
which is introduced into the definition of the eddy- viscosity. Both models are tested for a large number of 
different flowfields. The results of the BSL model are similar to those of the original k — u> model, but 
without the undesireable freestream dependecy. The predictions of the SST model are also independent of 
the freestream values and show excellent agreement with experimental data for adverse pressure gradient 
boundary-layer flows. 

INTRODUCTION 

The main field of application of Navier-Stokes methods in aerodynamics will be for complex turbulent 
flows that cannot be treated by inviscid, or viscous-inviscid interaction schemes. Examples are massively 
separated flows, flows involving multiple length-scales, flows with three-dimensional separation and 
complex unsteady flows. In these flows, the application of algebraic turbulence models, like the Cebeci- 
Smith [1], the Baldwin-Lomax [2] or the Johnson-King [3] model, becomes very complicated and often 
ambiguous, mainly because of the difficulty to define an algebraic length-scale. It is obvious that the 
improvement of numerical methods must be accompanied by the development of more general turbulence 
models and their implementation into Navier-Stokes codes. 

In addition to being independent of the specification of an algebraic length-scale, there is a long 
wish list of characteristics a good turbulence model would have to satisfy. Obviously, the model should 
be "sufficiently" accurate for the intended type of applications. Furthermore, and almost as important, 
the model must be numerically robust and should not consume excessive amounts of computation time 
(compared to the mean-flow solver). Another important demand is that the results should not have a strong 
dependency on ambiguous quantities, especially on the specified freestream values. 

The simplest nonalgebraic turbulence models available are the one-equation models developed re- 
cently by Baldwin and Barth [4] and by Spalart and Allmaras [5]. Both models solve one partial differential 
equation for the eddy- viscosity and appear to be very robust numerically (note however, that the numerical 
procedures recommended by the authors are only first-order accurate and an extension to second-order 
seems desirable). The present author has included the Baldwin-Barth model into a comparative study of 
Hiff prent eddy-viscosity models under adverse pressure-gradient conditions [6] and found that this model 
did predict the overall displacement effect of the boundary-layer (and its influence on the pressure dis- 
tribution) accurately, but failed to reproduce the proper shape of the velocity profiles. Furthermore, the 
model predicts significantly too low spreading rates for free shear layers, if small freestream values are 
specified for the eddy viscosity. The newer model of Spalart and Allmaras [5] predicts better results for 
free shear layers and more realistic velocity profiles. Unfortunately, it seems to be less sensitive to the 



effect of adverse pressure-gradients [5] and the present author has found that it gives very high backflow 
velocities in the recirculation region of a backward-facing step flow. Regarding the short time-span since 
the introduction of these models however, further improvements are certainly to be expected. 

The most popular nonalgebraic turbulence models are two-equation eddy-viscosity models. These 
models solve two transport equations, generally one for the turbulent kinetic energy and another one related 
to the turbulent length- (or time-) scale. Among the two-equation models, the k — t model is by far the 
most widely used today. The first low Reynolds number k - e model has been developed by Jones and 
Launder [7] and has subsequently been modified by many authors. 

The k - t model has been very successful in a large variety of different flow situations, but it also 
has a number of well known shortcomings. From the standpoint of aerodynamics, the most disturbing of 
them is the lack of sensitivity to adverse pressure-gradients. Under those conditions, the model predicts 
significantly too high shear-stress levels and thereby delays (or completely prevents) separation [8]. Rodi 
[9] attributes this shortcoming to the overprediction of the turbulent length-scale in the near wall region 
and has shown that a correction proposed by Hanjalic and Launder improves the predictions considerably. 
However, the correction is not coordinate-invariant and can therefore not be applied in general coordinates. 
An alternative way of improving the results has been proposed by Chen and Patel [10] and by Rodi [11]. 
They replace the e-equation in the near wall region by a relation that specifies the length-scale analytically. 
This also reduces some of the stiffness problems associated with the solution of the model. Although 
the procedure is coordinate independent, it has only been applied to relatively simple geometries, where 
the change between the algebraic relation and the e-equation could be performed along a pre-selected 
gridline. Clearly this cannot be done in flows around complex geometries. Furthermore, the switch has to 
be performed in the logarithmic part (the algebraic length-scale is not known in the wake region), so that 
the original k — t model is still being used over most of the boundary layer. 

Another problem with the k — t model is associated with the numerical stiffness of the equations when 
integrated through the viscous sublayer. This problem clearly depends on the specific version of the k - e 
model selected , but there are some general aspects to it. All low Reynolds number k - t models employ 
damping functions in one form or another in the sublayer. These are generally highly nonlinear functions 

of dimensionless groups of the dependent variables like Rt = (models involving ?/ + are undesirable in 

separated flows). The behavior of these functions cannot easily be controlled by conventional linearization 
techniques and can therefore interfere with the convergence properties of the scheme. A second problem 
is that e does not go to zero at a noslip surface. That in turn leaves two alternatives. One is to employ a 

nonlinear boundary condition on e (e = const. &)(£-), or to add additional terms to the e-equation [7] that 
allow the use of a homogeneous boundary condition. Both approaches introduce additional nonlinearities 
that can upset the numerical procedure. 

There is a significant number of alternative models [13, 14] that have been developed to overcome 
the shortcomings of the k - e model. One of the most successful, with respect to both, accuracy and 
robustness, is the k — u> model of Wilcox [12]. It solves one equation for the turbulent kinetic energy k 
and a second equation for the specific turbulent dissipation rate (or turbulence frequency) w. The model 
performs significantly better under adverse pressure-gradient conditions than the k - t model although 
it is the authors experience that an even higher sensitivity to strong adverse pressure-gradients would be 
desirable [6]. Another strong-point of the model is the simplicity of its formulation in the viscous sublayer. 
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The model does not employ damping functions and has straightforward Dirichlet boundary conditions. 
This leads to significant advantages in numerical stability. 

However, the k - u model also has a major shortcoming. It has been reported recently that the results 
of the model depend strongly on the freestream values, uj, that are specified outside the shear-layer. In [15] 
this problem has been investigated in detail, and it has been shown that the magnitude of the eddy- viscosity 
can be changed by more than 100% just by using different values for ujj. This is clearly unacceptable and 
corrections are necessary to ensure unambiguous solutions. Another point of concern is that the model 
predicts spreading rates that are too low for free shear-layers, if the correct values are specified for uj. 

In this paper, two new turbulence models will be presented. First, a new baseline (BSL) k — u> model 
will be described. It is identical to the k - u model of Wilcox [12] for the inner region of a boundary 
layer (up to approximately S/2) and gradually changes to the high Reynolds number version of the k — t 
model of Jones and Launder (JL) in the outer wake region. In order to be able to perform the computations 
with one set of equations, the Jones-Launder model was first transformed into a k — u formulation. The 
blending between the two regions is performed by a blending function F\ that gradually changes from 
one to zero in the desired region. No a priori knowledge of the flowfield is necessary to perform the 
blending. The function also ensures that the Jones-Launder formulation is selected for free shear-layers. 
The performance of the new (BSL) model is very similar to that of the Wilcox k — w model for adverse 
pressure gradient boundary-layer flows (and therefore significantly better than that of the k t model), but 
without the undesirable freestream dependency. For free shear layers the new model is basically identical 
to the Jones-Launder k — t model. 

Although the original - and the new BSL k-u model perform better in adverse pressure gradient 
flows than the it - t model, they still underpredict the amount of separation for severe adverse pressure 
gradient flows [6]. In an attempt to improve matters, the eddy-vicosity formulation of the BSL model will 
be modified to account for the transport effects of the principal turbulent shear-stress. The motivation for 
this modification comes from the Johnson-King (JK) model [3] which has proven to be highly successful 
for adverse pressure gradient flows. The JK-model is based on the assumption that the turbulent shear- 
stress is proportional to the turbulent kinetic energy in the logarithmic - and the wake region of a turbulent 
boundary layer. Johnson and King solve an equation for the maximum turbulent shear-stress at each 
downstream station and limit the eddy-viscosity in order to satisfy this proportionality. In the framework 
of two-equation models the turbulent kinetic energy is already known and it is therefore only necessary 
to limit the eddy-viscosity to account for the same effect. The resulting model will be called shear-stress 
transport (SST) model. First predictions based on this assumption have already been reported in [6]. 

The BSL and the SST models are not significandy more complicated than the original k-u model 
and consume only little more computing time. Furthermore, the models have shown the same numerical 
robustness as the original model for all the flows computed so far. 

THE TURBULENCE MODEL 
The new Baseline (BSL) Model 

The basic idea behind the BSL model is to retain the robust and accurate formulation of the Wilcox 
jfc - u> model in the near wall region, and to take advantage of the freestream independence of the k - e 
model in the outer part of the boundary-layer. In order to achieve this goal, the k - e model is transformed 
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into a k — u formulation. The difference between this formulation and the original k — u model is that 
an additional cross-diffusion term appears in the w-equation and that the modeling constants are different. 
The original model is multiplied by a function F\ and the transformed model by a function ( 1 - F\ ) and 
both are added together. The function F\ will be designed to be one in the near wall region (activating the 
original model) and zero away from the surface. The blending will take place in the wake region of the 
boundary-layer. 


Original k — u model: 
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dt dx j 
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Now, equation 1 and equation 2 are multiplied by F\ and equation 3 and equation 4 are multiplied by 
(1 — F\) and the corresponding equations of each set are added together to give the new model: 
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(7) 


( 8 ) 


Let <f > i represent any constant in the original model {<J k \, ...) , <j> 2 any constant in the transformed k - e 
model ( 0 * 2 , —) and <f> the corresponding constant of the new model ( cr k ...), then the relation between them 
is: 

<j> = F\<f>\ - F\)4>2- (9) 

the following two sets of constants will be used: 


Set 1 (Wilcox): 
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< 7*1 = 0.5, <t w i = 0.5, Pi = 0.0750, _ 

P* = 0.09, k = 0.41, 71 = PilP* 


Set 2 (Jones-Launder): 

< 7* 2 = 1.0, < 7 w2 = 0.856, p2 = 0.0828, __ 

/?* = 0.09, k = 0.41, 72 = fc/P* - ^ 2 « 2 /v^- 

Set 1 corresponds to the original k — u model and will be used in the near wall region exclusively. Set 2 
corresponds to the exact transformation of the Jones-Launder model (c\ e = 1 .44, C 2 e = 1 «92) and its main 
area of application is for free shear-layers. 

The model has to be supplemented by the definition of the eddy-viscosity: 


The turbulent stress tensor = — u(tt' is then given by: 


,du{ 9uj 2 jr 
Ti > = + W ~ 3 


In order to complete the derivation of the model it is necessary to define the blending function F { . 
Starting from the surface, the function should be equal to one over a large portion of the boundary layer 
in order to preserve the desirable features of the k — oj model, but go to zero at the boundary layer edge 
to ensure the freestream independence of the k — e model. The function will be defined in terms of the 


variable: 


, • , Vk 400 v, 

arg, = "PC') 


( 12 ) 


as follows: 

F\ = tanh(argi) (13) 

with n being the absolute value of the vorticity. The first argument in equation 12 is the turbulent length 
scale Lt = Vfc/(0.09w) (= k z / 2 /e), divided by the distance to the next surface, y. The ratio Lt/y is equal 
to 2.5 in the logarithmic region of the boundary-layer and goes to zero towards the boundary-layer edge. 
The second argument is an additional safe guard against the “degenerate” solution of the original k - ui 
model with small freestream values for u [15]. The third argument in equation 12 simply ensures that the 
function F\ does not go to zero in the viscous sublayer. The behavior of this function can be controlled 
by multiplying L t /y by a factor (in this case 1) and by the exponent in equation 13. Figure 1 shows the 
typical behaviour of the function F\ for different velocity profiles in a strong adverse pressure gradient 
boundary layer. Figure 1 also includes the underlying velocity profiles (same linestyle). The function is 
equal to one over about 50% of the boundary-layer and then gradually goes to zero. The behavior of the 
new BSL model will obviously lie somewhere between the original k - u and the k - e model. However, 
since most of the production of both, k and w, takes place in the inner 50% of the layer, it can be expected 
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that the model performance will be closer to that of the k - model, governing this area. Recall that the 
replacement of the t equation by an algebraic length-scale, as proposed by [10, 1 1] has to be performed in 
the logarithmic region so that the original k — t model still covers the largest part of the boundary layer 
and their results will therefore be much closer to those of the A; — e model. 



Figure 1. Blending function F\ versus y/6 for different velocity profiles. 


The Shear-Stress Transport (SST) Model 

One of the major differences between an eddy-viscosity - and a full Reynolds-stress model, with 
respect to aerodynamic applications, is that the lat ter does account for the important effect of the transport 
of the principal turbulent shear-stress r =: -u'v' (obvious notation) by the inclusion of the term 


Dt 

Dt 


dr 


dr 


dt +Uk ~dx k 


(14) 


The importance of this term has clearly been demonstrated by the success of the Johnson-King (JK) 
[3] model. Note that the main difference between the JK - model and the Cebeci-Smith model lies in 
the inclusion of this term in the former, leading to significantly improved results for adverse pressure 
gradient flows. The JK model features a transport equation for the turbulent shear-stress r that is based 
on Bradshaw’s assumption that the shear-stress in a boundary-layer is proportional to the turbulent kinetic 
energy, k, : 

t = pa\k (15) 

with the constant a\ = 0.3. On the other hand, in two-equation models, the shear-stress is computed from: 

r = fi t Sl (16) 

with fl = For conventional two-equation models this relation can be rewritten to give: 


__ / Production k 

y Dissipation^ 


(17) 
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as noted in [6]. In adverse pressure gradient flows the ratio of production to dissipation can be significantly 
larger than one, as found from the experimental data of Driver [16], and therefore equation 17 leads to an 
overprediction of t. In order to satisfy equation 15 within the framework of an eddy-viscosity model, the 
eddy-viscosity would have to be redefined in the following way: 



( 18 ) 


The rational behind this modification can also be explained in the following way: In conventional two- 
equation models the turbulent shear-stress responds instantaneously to changes in the shear-strain rate ft, 
much like an algebraic eddy- viscosity model, whereas equation 18 guarantees that t does not change more 
rapidly than pa\k. Obviously, equation 18 is not desirable for the complete flowfield, since it leads to 
infinitely high eddy-viscosities at points where ft goes to zero. Note however, that in adverse pressure 
gradient flows, production is larger than dissipation for the largest part of the layer (or ft > a\u). The 
following expression: 

vt = ^ (19) 

max(a\u>]H) 

guarantees therefore the selection of equation 18 for most of the adverse pressure gradient regions (wake 
region of the boundary layer), whereas the original expression equatio n 10 is used for the rest of the 
boundary layer. Figure 2 shows the relation of (-u'v'/aik) versus ^{Production/ Dissipation) for 
the SST model (equation 19), the conventional k — u> {k — e) model (equation 10), Bradshaw s relation 
(equation 15) and a relation proposed by Coakley [13]. Note that Coakley’s relation contains the relations 
of Bradshaw (/3 = 1) and that of the conventional two-equation models (/3 = 0) as a subset, but not 
equation 19 (/? is a free parameter of Coakley’s model). 



In order to recover the original formulation of the eddy-viscosity for free shear-layers (where Brad- 
shaws assumption equation 15 does not necessarily hold) the modification to the SST model has to be 
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limited to boundary-layer flows. This can be done in the same way as for the BSL model by applying a 
blending function F 2 . 

a ]k 

n = 7 ttttt (20) 

where F 2 is defined similarly to equation 13: 


ar<72 = max (2 


Vk 400 V 

0.09 u)y' xpu 


) 


( 21 ) 


F<i = tanh(argl ) (22) 

F<i is depicted in figure 3 in the same way as F\ in figure 1. Since the modification to the eddy- viscosity 
has its largest impact in the wake region of the boundary layer, it is imperative that F 2 extends further out 
into the boundary-layer than Fj. (Note on the other hand that F\ has to fall off to zero well within the 
boundary-layer in order to prevent the freestream dependence of the k — u model). 



Figure 3. Blending function F 2 versus y/8 For different velocity profiles. 


This modification to the eddy-viscosity is used in connection with the BSL model derived above. 
However, the SST model produces slightly lower eddy-viscosities than the BSL model for a flat plate zero 
pressure gradient boundary layer. In order to recover the correct cy-distribution the diffusion constants in 
the near wall model had to be adjusted accordingly: 

Set 1 (SST - inner): 


cr/ci = 0.85, <r u 1 = 0.65, = 0.0750, 

0 * = 0.09, k = 0.41, 71 = A//9* - 
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Set 2 remains unchanged. Furthermore, for general flows f l is taken to be the absolute value of the vorticity. 
Both models are given in their full form in Appendix A. 

Two new turbulence models have been introduced in this section. Both can be regarded as zonal 
models, since they utilize different models for different areas of the flowfield. However, in contrast to 
the classical zonal approach the present models do not require an a priori knowledge of the flowfield in 
order to define the zonal boundaries where the different models are to be used. The change between the 
different sub-models is achieved by "smart" functions that can distinguish between the different zones. 
Note that similar functions can also be designed for the k — t model. This makes it possible to design one 
set of constants for boundary-layer flows and a second set for free shear-layers and to switch between the 
different sets in the same way as in equation 9. 

The BSL model retains the advantages of the Wilcox k - to model for boundary-layer applications, 
but avoids its undesirable freestream dependency. Furthermore it switches to the more accurate k — t 
model for free shear-layer applications. In addition to this, the SST model modifies the definition of the 
eddy-viscosity for adverse pressure gradient boundary-layer flows in much the same way as the Johnson- 
King model does. From a computational point of view both models are not significantly more complex 
than the original k — to model. 

Boundary Conditions 

At a no-slip wall, all turbulent quantities, except to are set to zero. As pointed out by Wilcox [ 12], w 
satisfies the following equation near the wall: 


to 


6v 

P iP 


as 


0 . 


(23) 


Wilcox recommends to specify this analytical solution for the first few gridpoints away from the wall 
explicitly. The present author found it much easier and as accurate to implement the following boundary 
condition: 

" = 10 ^ « »=° <24) 


where A y is the distance to the next point away from the wall. Equation 24 simulates the boundary 
condition 23 without the need of changing the solution at interior points. It should be noted that models 
based on the ^-equation give accurate results if the near wall values of to are sufficiently large. Both, 
equations 23 and 24 satisfy this demand. The results are not sensitive to the factor ( 1 0) used in equation 24. 


At inflow boundaries, the turbulence quantities are specified and at outflow boundaries a zero gradient 
is assumed. 


Two of the computed flowfields have a rotational symmetry. In this case, the gradients of all quantities 
in the circumferential direction are set to zero. 

NUMERICAL METHOD 


The mean flow equations are solved by the INS3D code of S. E. Rogers and D. Kwak [18] which is 
based on a pseudo-compressibility method. The convective terms are upwind differenced with a third-order 
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flux-difference scheme. The viscous fluxes are differenced with second-order accurate central differences. 
The linear equations resulting from the first-order backward Euler time differencing are solved with an 
implicit line relaxation scheme. 

The turbulence equations are solved by a very similar method. However, instead of the third-order 
upwind differencing, the second-order TVD (total variation diminishing) scheme of [19] is applied. It 
was found that the TVD scheme prevents small oscillations of k and u> near the boundary-layer edge and 
thereby improves the stability of the scheme considerably. The turbulence equations are solved decoupled 
from the mean flow equations in a separate subroutine. One of the important aspects in the discretization of 
turbulence models is the implicit treatment of the source terms. The following approximate linearization 
gave good numerical properties: 

(25) 

^-{Pw +c u - D u ) * J£?) +2Dw . ( 26 ) 

ou> u > 

where P, D, C are the production, the destruction and the additional cross diffusion terms respectively. 
The above expressions go to the left hand side of the algorithme with a change of sign and thereby increase 
the diagonal dominance of the algorithm. The resulting numerical scheme has proven to be very robust 
and all of the computations with the k — u> models could be run with an infinite time step. An exception 
is the backward facing step flow where the time step had to be reduced for all models tested so far by the 
author. Furthermore, the computations could be started with very crude initial conditions (like freestream 
values). 


All computations, except the airfoil computations, have been performed on different grids, to ensure 
that the presented solutions are grid independent. The airfoil compuations were performed on a standard 
grid, kindly provided by S. E. Rogers [20]. 

RESULTS 

Flat Plate Boundary Layer 


In order to show the motivation for the derivation of the BSL model, flat plate zero pressure gradient 
boundary-layer computations with different freestream values for u> have been performed. It has been 
shown in [15] that the correct freestream values for u> outside the boundary-layer are: 


w / = 


4 

y/FU e 6*' 


(27) 


where u T = \Jt w I p is the friction velocity, defined in terms of the wall shear-stress t w , U e is the freestream 
velocity and 6* is the displacement thickness. For the first set of computations, the above value has been 
specified at the inflow boundary in the freestream for both, the original and the BSL k — u> model. Then, 
the above value was reduced by four orders of magnitude and the computations were repeated with both 
models. Note that the freestream value of k was also reduced in order to keep the freestream value of the 
eddy-viscosity constant (equal to the molecular viscosity). Figure 4 shows eddy-viscosity profiles for the 
original and the BSL k — u model. The eddy- viscosity of the original model changes by almost 100% due 
to the changes in w*, whereas the BSL model gives the same results for both cases. The strong sensitivity 
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of the original model to Uf is clearly unacceptable and can lead to a severe deterioration of the results for 
complex flows, as will be shown later. Results of the SST model were also found to be independent of u>f. 




Figure 4. Freestream-dependency of the eddy-viscosity for the original and the BSL k — u model 

In each of the following comparisons between the different models, u>f for the original k — u model 
was always chosen according to equation 27 whenever possible. In cases where the ffeestream values had 
to be chosen differently, it will be mentioned in the text. Figure 5 shows a comparison of the SST, the 
BSL, the original it - u and the Jones-Launder (JL) k - e model (all JL model computations have been 
performed with damping functions as given in [?]) for a zero pressure gradient flat plate boundary layer. 
Obviously, all models predict the correct cy-distributions and velocity profiles. 

Computations have been performed on a grid with 90x90 gridpoints. 

Free Shear Layers 

For free shear layers the SST and the BSL model reduce to the same model (F 2 = 0; = 0), which 

will be called new k - u> model in this subsection. Note that the new k - u> model is formally almost 
identical to the Jones-Launder k - e model. However a small cross-diffusion term has been neglected in 
the derivation of equation 4. In order to show that this term is truly negligible, the original k - e model is 
also included in the comparison. All computations were performed with 200 gridpoints across the layer. 
A grid refinement study on a grid with 300 points gave the same results. 

Figure 6 shows the results of the solution of the equilibrium far wake equations for the different 
models. The results are compared to the experimental data of Fage and Falkner [24]. Obviously, the new 
k - u> model and the JL k - e model produce almost identical results. The original k - u model predicts 
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Re,/ 10" y + 


(a) Skin friction (b) Law-of-the-wall 

Figure 5. Flat plate boundary-layer. 

a somewhat lower spreading-rate than the other models. As is well known for two-equation models, they 
all fail to predict the smooth behavior of the data at the edge of the layer. The freestream value for the 
original k — u model has been derived from an expression similar to equation 27 [15]. 


1.0 
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0.6 
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0.2 
0.0 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 

Figure 6. Comparison of the velocity distribution of three different two-equation models for the far wake 
flow of Fage and Falkner 

Figure 7 shows solutions for a self-similar plane jet, as reported in [22]. Again, the new k — u and 
the JL model produce almost identical results in very good agreement with the experiments, whereas the 
original k — ui model shows a rather peculiar behavior. The major difference to the far- wake computations 
is that for the present flow the freestream velocity is zero (still air). The only acceptable freestream value 
for u> in still air is u>f = 0, as can be seen from the equilibrium equations [15] and from physical intuition. 
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The specification of small values for uj leads to large eddy- viscosities in the original k - w model, as 
demonstrated in figure 4 for the flat plate boundary-layer. The same is true in the present case to an even 
larger extent, because of the missing wall influence. The original k — u> model predicts about five times as 
high an eddy-viscosity as the other two models, resulting in the large spreading rates shown in figure 7. 



Figure 7. Comparison of the velocity distribution of three different two-equation models for the plane jet 
flow of Wygnanski and Fiedler 


The last free shear-layer in this comparison is the free mixing-layer of Liepman and Laufer [23], 
figure 8. Note that the freestream velocity below the mixing-layer is zero, leading to the same problem 
with the original k — ui model as for the plane jet. The other two models again produce almost identical 
results in acceptable agreement with the experiments. 



Figure 8. Comparison of the velocity distribution of three different two-equation models for the free 
mixing-layer of Liepman and Laufer 
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Adverse Pressure Gradient Flows 


One of the most important aspects of a turbulence model for aerodynamic applications is its ability to 
accurately predict adverse pressure gradient boundary-layer flows. It is especially important that a model 
can predict the location of flow separation and the displacement effect associated with it. The reason is that 
the viscous-inviscid interaction has a strong influence on the overall pressure distribution and therefore on 
the performance of the aerodynamic body. 

The most widely used test case to measure the performance of turbulence models under adverse 
pressure gradient conditions is the flow reported by Samuel and Joubert [25]. It is a flat plate boundary- 
layer, developing under an increasingly adverse pressure gradient The upstream Reynolds number is 
1 .7 • 10 6 m -1 . The flow is retarded by the pressure gradient, but not separated. 

Two different sets of computations have been performed. At first, the experimental pressure distribu- 
tion was specified at the outer edge of the computational domain (opposite to the wall). Since this is not 
a very straightforward boundary condition for a Navier-Stokes method, a second set of computations was 
performed based on the specification of an inviscid external streamline. The inviscid streamline y(x) s is 
defined by the preservation of mass: 


y(x) s 

m = J pu eX p.{x,y)dy = const. (28) 

0 

Note that the the specification of a streamline does not mean that the displacement thickness is prescribed 
like in an inverse boundary-layer method. Both computations produced very similar results. The solutions 
shown here are for the prescribed streamline, which is thought to be the more consistent boundary 
condition from a numerical point of view. The eddy-viscosity at the inflow boundary was determined from 
the exp erim ental shear-stress and velocity profiles, the turbulent kinetic energy k from the requirement 
k = (—u'v')/ai and u> and e from the definition of the eddy-viscosity. The same grid of 90x90 points as 
in [6] was used for the compuations. The results are virtually identical to those on a 60x60 grid. 

Figure 9 shows a comparison of the computed and the experimental skin-friction distribution. All 
three k — u models reproduce the experimental data almost exactly, whereas the JL k — t model gives 
significantly too high values. 

Figure 10 shows the same comparison for the velocity profiles. There are no large differences to be 
found between the different models. Only the SST model predicts a somewhat stronger retardation of the 
flow near the surface. The same behavior has been observed with the Johnson-King model in [6] 

Turbulent shear-stress profiles are shown in figure 11. All models are slightly overpredicting the 
amount of shear-stress, with the SST model closest and the JL model furthest away from the data. 

It is obvious that the small differences between the solutions, especially between the different k - u 
models, does not allow final conclusions about the abilities of the models to predict adverse pressure 
gradient flows. It appears that the Samuel- Joubert flow does not pose a sufficiently strong challenge to 
the models to assess their potentials for this type of flows. The author has reached a similar conclusion 
in [6], where a solution of the Johnson-King (JK) model has shown, that the model did not significantly 
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Figure 11. Turbulent shear-stress profiles for Samuel-Joubert flow at x= 1.16, 1.76, 2.26, 2.87, 3.40 m. 

depart from its equilibrium formulation. It is therefore important to test models under more demanding 
conditions, with stronger adverse pressure gradients and possibly separation. The following flowfield, 
reported by Driver [16], has proven to be a very self-consistent and demanding test case, and is therefore 
strongly recommended for the assessment of turbulence models under adverse pressure gradients. 

In Driver’s flow, a turbulent boundary-layer develops in the axial direction of a circular cylinder. An 
adverse pressure gradient is imposed by diverging wind tunnel walls and suction applied at these walls. 
The pressure gradient is strong enough to cause the flowfield to separate. The inflow Reynolds number is 
2.8 • 10 5 based on the diameter D of the cylinder (140mm). 

The boundary conditions for this flow are similar to those used for the Samuel-Joubert flow. Again 
an inviscid streamline is extracted from the experimental velocity profiles and a slip condition is applied 
along this line. The inflow conditions are determined from the experimental profiles in the same way as 
described above. The computations have been performed with a three-dimensional version of the code. 
In order to account for the axial symmetry, three closely spaced circumferential planes where introduced 
and symmetry conditions were applied. A 60x3x60 grid [6] was used for the present computations. A 
computation on a 100x3x100 grid gave almost identical results. 

Figure 12 shows the wall pressure distribution for Driver’s flow as computed by the different models. 
The SST model gives superior results to the other models due to its ability to account for the transport of 
the principal shear-stress. As expected, the JL k - t model produces the worst results, and the BSL and 
the original k — u model being close to each other in the middle. 
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Figure 12. Wall pressure distribution for Driver’s adverse pressure-gradient flow. 

Figure 13, depicting the wall shear-stress distribution for Driver’s flow, shows that the SST model 
predicts the largest amount of separation, whereas the JL model stays firmly attached. Again, the BSL and 
the orig. k — u> model produce very similar results, in good agreement with the experiments. 



x/D 


Figure 13. Wall shear-stress distribution for Driver’s adverse pressure-gradient flow. 


The differences between the models can be seen more clearly in figure 14 which shows the velocity 
profiles. The SST model clearly produces the best agreement with the experiments. The larger displacement 
effect predicted by this model is reflected in the flattening of the c p -distribution as observed in figure 12. 
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The orig. k — u> model predicts slightly better results than the BSL model, and the JL k — t model shows 
very little sensitivity to the pressure gradient, as already reflected in figures 12 and 13. 



Figure 14. Velocity profiles for Driver’s adverse pressure-gradient flow at x/D=-0.09 1,0.363, 1.088, 1.633, 
2.177. 


The reasons for the different behavior of the models can be seen in the following two pictures. 
Figure 15 compares turbulent shear-stress profiles at different stations. The experimental data are shown 
both, in a surface (Carth.) oriented and in a streamline (Strml.) oriented coordinate system. (Note that the 
numerical results were, due to the eddy-viscosity formulation, not found to be sensitive to small changes 
in direction). Driver [17] has shown that the streamline direction is perpendicular to the direction of the 
maximum shear-strain rate for this flow. Since eddy-viscosity models are calibrated in this direction it is 
the authors opinion that this is the most meaningful direction for the comparison. 

The JL model obviously predicts significantly higher shear-stress levels than the other models, 
especially in the region where separation is approached. This in turn leads to the firmly attached velocity 
profiles of figure 14. The differences between the models can be seen more clearly by looking at the 
eddy-viscosity distributions. Figure 16 shows the maximum value of the kinematic eddy-viscosity profiles 
for all x-stations, nondimensionalized by u e S*. The SST model predicts the reduction of this quantity due 
to the adverse pressure gradient in very good agreement with the experiments. The BSL and the orig. k-u 
model are very close to each other up to separation (around x/D = 0), whereas the orig. model is closer to 
the experiments in the recovery region. Both these models are giving values that are consistently too high. 
The k - t model falls only barely below the value of 0.0168 recommended by Clauser for equilibrium 
b idary-layers (and used in the Cebeci-Smith model) and thereby fails to account for the nonequilibrium 
enacts altogether. 
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Figure 15. Turbulent shear-stress profiles for Driver’s adverse pressure-gradient flow at x/D=-0.09 1 , 0.363, 
1.088, 1.633, 2.177. 

Backward-Facing Step Flow 

Results for the flow over a backward facing step as reported by Driver and Seegmiller [26] will 
be discussed next. This flowfield was a test case in the 1981 Stanford conference for the evaluation of 
turbulence models. However, most of the computations at the time have been performed on comparatively 
coarse grids and there is substantial evidence that significantly finer grids have to be used in order to obtain 
grid-independent results [27]. The present computations have been performed on a 120x120 grid, with 
substantial grid refinement near the step. Figure 17 shows the distribution of gridpoints. As with the other 
flowfields (except the airfoil, where a standard grid was used), a grid refinement study was made. The 
present results are virtually identical to those performed on a 90x90 grid. The Reynolds number, based on 
the upstream momentum thickness 0 is, = 5000 and the ratio of the boundary-layer thickness to the 
step height is about 1.5. The expansion ratio (height of the tunnel after the step divided by the height of 
the tunnel in front of the step) is 1.125. 


Figure 18 shows a comparison of computed and experimental skin friction distributions. The k - u 
models all perform significantly better than the k - e model. The reattachment length of the four models are 
6.7 (SST), 5.9 (BSL), 6.4 (org. ku) and 5.5 (JL k - e) compared to a value of about 6.4 in the experiments. 
The original k - u model reproduces the reattachment location almost exactly, whereas the SST and the 
BSL model are very close to the experiments. The reattachment length predicted by the k - e model is 
better than previously reported, certainly as a result of the fine grid employed in the present computations 
(see also [27]). However, the model predicts significantly too large variations of cj in the recirculation and 
the reattachment region. The good results of the fc - w models show that it is not necessary to account for 
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the anisotropy of the stress tensor, as suggested in [28], in order to get accurate results for the reattachment 
length. 



Figure 18. Wall shear-stress distribution for backward-facing step flow. 


The surface pressure distribution, shown in figure 19, reflects the trends already seen in cj. The larger 
the separation region predicted by the model the larger is the displacement effect of the boundary layer and 
the smaller is the pressure rise in the expansion region after the step. The SST model slightly overpredicts 
separation and therefore gives a somewhat smaller pressure increase than the other models. 

Figure 20 shows a comparison of the velocity profiles at x/H=2.0, 4.0, 6.5, 8.0, 14.0, 32.0. The figure 
shows again that the SST model predicts a larger separation region than the other models. All models fail 
to capture the relaxation downstream of reattachment correctly. The results of [28] show that this is also 
true for more complex models which account for anisotropy effects. 

Although the SST and the BSL models give slightly less accurate results than the orig. k - u> model, 
the overall agreement with the experimental data is still remarkable, considering the complexity of the 
problem. 

NAGA 4412 Airfoil Flow 

As an external flow test of the turbulence model, computations were performed for the flow around 
a NACA 4412 airfoil at 13.87 degrees angle of attack. The Reynolds number with respect to the chord 
length is Re = 1 .52 • 10 6 . Experimental data for this flow have been reported by Coles and Wadcock [29]. 
The grid for the computations consists of 241x61 gridpoints and was made available by S. Rogers [20]. It 
is similar to the one used in [30]. 

In the experiment the transition was fixed by transition strips at x/c of 0.023 and 0.1 for the upper 
and lower surface respectively. As reported in [30], if these locations are specified in the computations, 


21 




Figure 19. Wall pressure distribution for backward-facing step flow. 



Figure 20. Velocity profiles for backward-facing step flow at the streamwise locations: x/H= 2.0, 4.0, 6.5, 
8.0, 14.0, 32.0. 
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a laminar separation bubble appears before the transition point on the upper surface of the airfoil. This 
separation bubble was not observed in the experiments which indicates that transition may take place 
already before the strip is reached. Computations have been performed with and without a specified 
transition location and differences between the computations are small. Results are given here for the case 
were transition was not specified, so that the models picked their own transition location. 

In the computations without a specified transition location an important observation was made in 
the vicinity of the stagnation point. All three k — u models produced very high eddy-viscosities in this 
area, even outside the boundary layer. This in turn produced results far from the experimental data. The 
reason for this unphysical behavior comes from the incorrect prediction of the normal stresses by the 
eddy- viscosity model. The production of the tubulent kinetic energy k in regions of zero shear-strain (like 
in the stagnation area) is proportional to: 


P k — const.(—u ,2 + 


(29) 


In an eddy-viscosity model this translates to: 


P k = 4 const.nt{—) 2 . 


(30) 


Note that the latter expression is positive for accelerated and decelerated flows alike, which is clearly 
unphysical. The failure of the eddy-viscosity models in regions of dominating normal strain rates is not 
surprising, since these models are calibrated with respect to the shear strain. In order to prevent this 

undesireable behavior, the group (f^Xf^ + in the production of k and u is replaced by the square 
of the vorticity Cl: 

P(P + (31) 

OX j OXj C/Xj 

The effect of this modification is that the production of eddy-viscosity in inviscid flow regions is avoided. 
After this modification transition takes place inside the boundary-layer at a downstream station of x/c « 
0.006 on the upper and x/c « 0.06 on the lower surface of the airfoil. The modification has also been 
tested for the other flows presented above and produces almost identical results (even for the backward 
facing step). 

Figure 21 shows a comparison of the computed and the experimental velocity profiles at different 
stream wise stations. The results are similar to those for the separated case of Driver, figure 14. Again, 
the SST model predicts the displacement effect in very good agreement with the experiments. The BSL 
model is showing some response to the pressure gradient, and produces results similar to those reported in 
[30] for the Baldwin-Barth model. Another interesting result of this computation is that the original k - u> 
model predicts velocity profiles even further away from the experiments than the Jones-Launder k — e 
model. The reason for the poor performance of the orig. k - u model lies in its freestream dependency. 
In order to understand the problem, one has to look at the development of u> from the inflow boundary to 
the leading edge of the airfoil. In this inviscid region production and diffusion of u are zero. The balance 
in the u> equation reduces therefore to: 


U 8 ^ = -/3a; 2 
as 


(32) 
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Figure 21. Velocity profiles on the upper surface of a NACA 4412 airfoil at 13.87 degrees angle of attack. 
Streamwise stations x/c= 0.675, 0.73 1 , 0.786, 0.842, 0.897, 0.953. 


where s is the streamline direction and U s is the velocity in this direction. Assuming that (J s is constant 
and equal to Uqo, the equation can be solved to give: 


w(s) 



(33) 


The largest value of u that can be achieved for a certain distance s from the inflow boundary (s = 0) is: 


w(s) = 



(34) 


corresponding to an infinitely large u>oo- As s gets larger this maximum value becomes smaller and smaller. 
In the present computations the distance between the inflow boundary and the airfoil is about 15 chord 
lengths. Non-dimensionalizing all quantities with Uqo and the chord length, c, leads to a freestream value 
of u> in the leading edge region of the airfoil of about uj = 1 whereas the formula given in [15] for the 
estimation of the correct freestream value: 


Wf = 


ut 


VPU e 6* 


(35) 


indicate^ at it should be about three orders of magnitude larger. The low freestream value of a; in turn 
leads to me very high eddy-viscosities shown in figure 4 which in turn prevent separation. This example 
clearly shows the dangers of using the orig. k - u model for aerodynamic applications. If the correct 
freestream values could be specified, the results of the orig. k-u> should be very close to that of the BSL 
model. 


Figure 22 shows a comparison of the computed and the experimental surface pressure distributions. 
The agreement, especially for the SST model is not as good as has to be expected from the velocity profiles 
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shown in figure 21. Although the SST model predicts the displacement effect of the boundary layer 
almost exactly, it fails to reproduce the pressure distribution. This indicates an inconsistency between the 
computations and the experiments. Likely candidates are blockage effects in the wind tunnel (wind tunnel 
walls were not included in the computations) or three-dimensional effects in the experiment. 



Figure 22. Surface pressure distribution for a NACA 4412 airfoil at 13.87 degrees angle of attack. 

CONCLUSIONS 


Two new turbulence models have been derived from the original k-u model of Wilcox [12]. The 
motivation behind the new baseline (BSL) model was to eliminate the freestream dependency of the k - u 
model but retain its simple and reliable form, especially in the near wall region. In order to achieve this 
goal, a switching function was designed that can discriminate between the inner part (appr. 50%) of a 
boundary-layer and the rest of the flowfield. In this inner part the original k-u model is solved, and in 
the outer wake region a gradual switch to the high Reynolds number version of the Jones-Launder k - t 
model, in a k — u formulation, is performed. 

The BSL model was then used to derive a model that can account for the transport of the turbulent 
shear stress (Shear-Stress Transport or SST model). The derivation of the model was inspired by the 
success of the Johnson-King (JK) model. The equilibrium model underlying the JK formulation performs 
similar to a standard k - t model for adverse pressure gradient flows. By including the effect of the 
“inertia” of the transport terms into a transport equation for the principal shear-stress the model is able 
to produce highly accurate results for a large variety of flow problems. (Note that the improvement for 
adverse pressure gradient flows comes from a change in the outer (wake-region) part of the eddy-viscosity 
and not from a more refined modeling in the sublayer). The main assumption in the JK model that the 
principal shear-stress is proportional to the turbulent kinetic energy was incorporated into the new SST 
model. This modification ensures that the principal shear-stress satisfies the same transport equation as the 
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turbulent kinetic energy. It is designed to act only inside the boundary-layer in order to retain the k — e 
model (transformed to a k — u> formulation) for free shear-layers. 

Both models were applied to a selection of well documented research flows, that are meaningful 
for aerodynamic applications. The results of the computations were compared against solutions of the 
standard k — u> and the standard k — t model, as well as against experimental data. 

The free shear-layer results have shown that the new models give results almost identical to those 
of the k — c model. Another important result of those computations is that they show clearly the strong 
ambiguity in the results of the original k — u model with respect to freestream values. 

The central part of the comparisons is for the behavior of the models under adverse pressure gradient 
conditions. The computations of the Samuel- Joubert flow, as well as Driver’s separated adverse pressure 
gradient flow show that the SST model gives highly accurate results for this type of problem. The BSL 
and the original k - u> model produce rather similar results, provided the correct freestream values are 
specified in the latter. 

Computations were also performed for the backward facing step flow of Driver and Seegmiller 
[26]. A very fine grid was employed to ensure grid independence of the results. For this problem, the 
original k - u> model gives very accurate results. It predicts the reattachment length within the uncertainty 
of the measurements and gives an accurate representation of the wall pressure distribution. The SST 
and the BSL models are giving about 5% too large or 8% too small values for the reattachment length 
respectively. These results are still very accurate, considering the notorious difficulties this flow poses to 
numerical assessment. All models fail to predict the relaxation of the velocity profiles downstream of the 
reattachment point correctly. The computations also show that it might not be necessary to account for 
anisotropy effects in the stress tensor, as suggested in [28], in order to predict the reattachment location 
correctly. The last set of computations has been performed for a NACA 4412 airfoil at an angle of attack 
near maximum lift condition. The computations confirm the findings of the adverse pressure gradient 
computations. The SST model predicts highly accurate velocity profiles, almost identical to those of the 
experiments. The BSL model has a smaller sensitivity to the adverse pressure gradient than the SST model 
and therefore predicts less retarded profiles. A very surprising result of the computations is that the original 
k - u model gives even less accurate solutions than the Jones-Launder k - t model. The reason for the 
failure of the model is again its freestream dependency. This computation clearly shows that the original 
k — u model cannot be applied unambiguously for industrial types of applications. 

The SST model in its present formulation is somewhat more sensitive to adverse pressure gradient 
conditions than the Johnson-King model when applied to the same flows [6], Since the Johnson-King 
model is known to give very accurate results for airfoil and wing flows, it might become necessary to 
fine-tune the SST model with reagard to those problems. This can easily be done by increasing the constant 
a\ in equation 20 slightly, in order to increase the eddy-viscosity for adverse pressure gradient flows. By 
increasing a\ the SST model moves more towards the BSL model and therefore predicts less separation. 
Note that the diffusion constants and have to be ajusted towards their BSL values as well, in order 
to recover the correct skin-friction for a flat plate zero pressure gradient boundary layer. 
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APPENDIX 

THE BASELINE (BSL) MODEL 
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For airfoil computations the following modification to the production terms is recommended (can also be 
used for other applications): 


dui . dui du j N 

— -( — - h -) 

dxj v dxj dx{ ’ 


fr (A-5) 

L /Jyj UJL j 

where ft is the absolute value of the vorticity. 

The constants <j> of the new model are calculated from the constants, fa, fa, as follows: 

4> = F\fa + (\ - F\)fa. (A-6) 

The constants of Set 1 (fa) are (Wilcox): 

<Jk\ = 0.5, <r w i = 0.5, 0\ = 0.0750, 

0* = 0.09, k = 0.41, 71 = 0i/0* ~ 


The constants of Set 2 (fa) are (Jones-Launder): 


<t *2 = L0, <t w2 = 0.856, 02 = 0.0828, 

0* = 0.09, k = 0.41, 72 = 0 2 /0 * - a^/VF- 

The function F\ is defined as follows: 

F\ = tanh(argf) 

with: 

z • f \fk u 400i/ 

0 m = °' 45 n ); * 

where y is the distance to the next surface. 

The eddy-viscosity is defined as: 

k 

vt = ~ 

U> 


(A-7) 

(A-8) 


(A-9) 
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THE SHEAR-STRESS TRANSPORT (SST) MODEL 


The SST model is identical to the above formulation, except that the constants, <j>\, have to be changed 
to: 

Set 1 (SST - inner): 

a kl = 0.85, <r w i = 0.65, /?i = 0.0750, 

(3* = 0.09, k = 0.41, 71 = 



(A- 10) 

(A- 1 1) 
(A- 12) 


Important detail !: 


In applying this model, it is important that the reader is aware of the following ambiguity in the 
formulation of the production term of u for the SST model. The definition of the production term of u is 


sometimes written as: 


Pu> = 


u> du{ 

1 k Tij dT j 


(A- 13) 


which introduces the nondimensional group vtj in front of the strain rate tensor. In the original and in the 
BSL model this group is equal to one and the two formulations for P u are therefore identical. This is not 
the case for the SST model because of equation A-10. The SST model has been calibrated with respect to 
equation A-4 and equation A- 13 should therefore not be used. 
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